    
/*********************************************************************
*                                                                    *
* F -- evaluates the function f(y,t), lots of physical imput on it   *
*                                                                    *
* Parameters:                                                        *
*       fields -- pointer to field_vector from where to extract (y,t)*
*       derivs -- pointer to field_vector where to put derivatives   *
*     wave_par -- pointer to parameter struct                        *
*                                                                    *
* Returns: nothing                                                   *
*                                                                    *
*********************************************************************/
#include "first_macro_2d.h"  /* Where global parameters are defined */
#include "structs_2d.h"      /* Where structures are defined */
#include "derivs_2d.h"       /* Where derivatives functions are defined */
 
#include "gen_2d.h"

#define GRD_IND [IND_2(grid_ind1,grid_ind2)]


void wave_2d_eq(struct GRID_PAR *grid_ptr,
		struct field_array *fields_ptr, 
		struct field_array *derivs_ptr,
		struct FUNCTION_PAR *function_par) {


    int ni_1 = (*grid_ptr).start_grid_1; 
    int nf_1 = (*grid_ptr).final_grid_1;  
    int ni_2 = (*grid_ptr).start_grid_2;  
    int nf_2 = (*grid_ptr).final_grid_2; 
    int n_gridpts_1 = nf_1 - ni_1;
    int n_gridpts_2 = nf_2 - ni_2;
    int grid = (*grid_ptr).grid;
    int coorx = 1;    
    int coory = 2;


    FLOAT xi = (*grid_ptr).initial_x;
    FLOAT xf = (*grid_ptr).final_x;
    FLOAT yi = (*grid_ptr).initial_y;
    FLOAT yf = (*grid_ptr).final_y;


    FLOAT one_dN1 = (xf-xi)/(FLOAT)(nf_1-ni_1-1);
    FLOAT one_dN2 = (yf-yi)/(FLOAT)(nf_2-ni_2-1);

    FLOAT sigma = (*function_par).sigma; /* For dissipation */

    FLOAT mm = (*function_par).mm;
    FLOAT ee = (*function_par).ee; 
    FLOAT d2, d2_aa, d2_bb, d2_cc, d2_ab, d2_cb, d2_ac, dip, dip_t ;
 
    /* normals */

    FLOAT nx_x0 = (*grid_ptr).nx[0][0];
    FLOAT ny_x0 = (*grid_ptr).ny[0][0];
    
    FLOAT nx_x1 = (*grid_ptr).nx[0][1];
    FLOAT ny_x1 = (*grid_ptr).ny[0][1];
    
    FLOAT nx_y0 = (*grid_ptr).nx[1][0];
    FLOAT ny_y0 = (*grid_ptr).ny[1][0];
  
    FLOAT nx_y1 = (*grid_ptr).nx[1][1];
    FLOAT ny_y1 = (*grid_ptr).ny[1][1];
    FLOAT nx;
    FLOAT ny;

    FLOAT m_a=0.0;
    FLOAT m_a_t=0.0;
    FLOAT m_a_x=0.0;
    FLOAT m_a_y=0.0;

    FLOAT m_b=0.0;
    FLOAT m_b_t=0.0;
    FLOAT m_b_x=0.0;
    FLOAT m_b_y=0.0;

    FLOAT m_c=0.0;
    FLOAT m_c_t=0.0;
    FLOAT m_c_x=0.0;
    FLOAT m_c_y=0.0;

    FLOAT norm_m=0.0; 
	FLOAT m_t_t_dot_m, L_x, L_y, L_t, L_up_x, L_up_y, m_t_dot_m ;

    FLOAT m_a_t_n=0.0;
    FLOAT m_a_x_n=0.0;
    FLOAT m_a_y_n=0.0;
    
    FLOAT m_b_t_n=0.0;
    FLOAT m_b_x_n=0.0;
    FLOAT m_b_y_n=0.0;

    FLOAT m_c_t_n=0.0;
    FLOAT m_c_x_n=0.0;
    FLOAT m_c_y_n=0.0;

    FLOAT m_a_t_t=0.0;
    FLOAT m_a_x_t=0.0;
    FLOAT m_a_y_t=0.0;

    FLOAT m_b_t_t=0.0;
    FLOAT m_b_x_t=0.0;
    FLOAT m_b_y_t=0.0;

    FLOAT m_c_t_t=0.0;
    FLOAT m_c_x_t=0.0;
    FLOAT m_c_y_t=0.0;

    FLOAT Sdet, g_up_xx, g_up_xy, g_up_yy, n_up_x, n_up_y, norm_n, norm_m_2;
 
    FLOAT SP, nV_a, nV_b, nV_c, V_a, V_b, V_c, Tau_p, V_dot_m;

    FLOAT  factor;

    char macro_value_strg[100];


     
    GET_MACRO_VALUE(DERIV);
    if (strcmp(macro_value_strg,"derivS_2d")==0) {
      factor = 2.0;
    }
    else if (strcmp(macro_value_strg,"derivQ_2d")==0) {
      factor = 48.0/17.0;
    }
    else if (strcmp(macro_value_strg,"derivQ_3_2d")==0) {
      factor = 11.0/3.0;
    }
    else if (strcmp(macro_value_strg,"deriv_strand_third_order_boundaries_sixth_interior_2d")==0) {
      factor = 43200.0/13649.0;
    }
    else if (strcmp(macro_value_strg,"deriv_strand_fourth_order_boundaries_eight_interior_2d")==0) {
      factor = 5080320.0/1498139.0;
    }
    else {
      factor=2.0; printf("check factor por penalty!!!%s!!!",macro_value_strg);
    }


 
    Tau_p = factor*(FLOAT)(nf_1 - ni_1 -1)/(xf-xi);

    FLOAT A = 1.0;

    /* ---------------------------------- PMI part --------------------------------------------- */

#include"mpi_SR.c"

 /*   --------------------------------------------------------------------------------------------- */


  /* first the time */

  (*derivs_ptr).time = (*fields_ptr).time;



  /* MAKE swap fields AND take derivatives */

#ifdef DEBUG
 printf("in eqn. before taking derivatives \n");
#endif




#ifdef FORCED_NORMALIZATION
// forced normalization, does not seems to work ok....
  {register int grid_ind1, grid_ind2;
  for (grid_ind1 = ni_1; grid_ind1 < nf_1; ++grid_ind1){
    for (grid_ind2 = ni_2; grid_ind2 < nf_2; ++grid_ind2){
 
        	m_a   = (*fields_ptr).u[M_A  ]GRD_IND;
        	m_b   = (*fields_ptr).u[M_B  ]GRD_IND;
        	m_c   = (*fields_ptr).u[M_C  ]GRD_IND;
 
 	norm_m = sqrt(m_a*m_a + m_b*m_b + m_c*m_c);

         m_a = m_a/norm_m;
         m_b = m_b/norm_m;
         m_c = m_c/norm_m;
 
 	(*fields_ptr).u[M_A]GRD_IND = m_a;
 	(*fields_ptr).u[M_B]GRD_IND = m_b;
 	(*fields_ptr).u[M_C]GRD_IND = m_c;
 }}}
#endif



 {register int grid_ind1, grid_ind2;
 for (grid_ind1 = ni_1; grid_ind1 < nf_1; ++grid_ind1){
   for (grid_ind2 = ni_2; grid_ind2 < nf_2; ++grid_ind2){
 
     	g_up_xx = auxfields.u_aux[G_UP_XX]GRD_IND;
	g_up_xy = auxfields.u_aux[G_UP_XY]GRD_IND;
	g_up_yy = auxfields.u_aux[G_UP_YY]GRD_IND;
 
	  Sdet = 1./sqrt(g_up_xx*g_up_yy - g_up_xy*g_up_xy);

	 auxfields.u_aux[SWAP_X]GRD_IND = 
	 (g_up_xx*(*fields_ptr).u[M_A_X]GRD_IND + g_up_xy*(*fields_ptr).u[M_A_Y]GRD_IND)*Sdet;
	 auxfields.u_aux[SWAP_Y]GRD_IND = 
	 (g_up_yy*(*fields_ptr).u[M_A_Y]GRD_IND + g_up_xy*(*fields_ptr).u[M_A_X]GRD_IND)*Sdet;
     }
   }
 }

DERIV(grid_ptr,(struct field *)&auxfields.u_aux[SWAP_X],
         (struct field *)&dfields.du[D_X_M_A_X],&coorx);
DERIV(grid_ptr,(struct field *)&auxfields.u_aux[SWAP_Y],
         (struct field *)&dfields.du[D_Y_M_A_Y],&coory);

{register int grid_ind1, grid_ind2;
 for (grid_ind1 = ni_1; grid_ind1 < nf_1; ++grid_ind1){
   for (grid_ind2 = ni_2; grid_ind2 < nf_2; ++grid_ind2){
 
     	g_up_xx = auxfields.u_aux[G_UP_XX]GRD_IND;
		g_up_xy = auxfields.u_aux[G_UP_XY]GRD_IND;
		g_up_yy = auxfields.u_aux[G_UP_YY]GRD_IND;

		Sdet = 1./sqrt(g_up_xx*g_up_yy - g_up_xy*g_up_xy);

		auxfields.u_aux[SWAP_X]GRD_IND = 
		(g_up_xx*(*fields_ptr).u[M_B_X]GRD_IND + g_up_xy*(*fields_ptr).u[M_B_Y]GRD_IND)*Sdet;
		auxfields.u_aux[SWAP_Y]GRD_IND = 
		(g_up_yy*(*fields_ptr).u[M_B_Y]GRD_IND + g_up_xy*(*fields_ptr).u[M_B_X]GRD_IND)*Sdet;
     }
   }
 }

DERIV(grid_ptr,(struct field *)&auxfields.u_aux[SWAP_X],
         (struct field *)&dfields.du[D_X_M_B_X],&coorx);
DERIV(grid_ptr,(struct field *)&auxfields.u_aux[SWAP_Y],
         (struct field *)&dfields.du[D_Y_M_B_Y],&coory);

{register int grid_ind1, grid_ind2;
 for (grid_ind1 = ni_1; grid_ind1 < nf_1; ++grid_ind1){
   for (grid_ind2 = ni_2; grid_ind2 < nf_2; ++grid_ind2){
 
		g_up_xx = auxfields.u_aux[G_UP_XX]GRD_IND;
		g_up_xy = auxfields.u_aux[G_UP_XY]GRD_IND;
		g_up_yy = auxfields.u_aux[G_UP_YY]GRD_IND;

		Sdet = 1./sqrt(g_up_xx*g_up_yy - g_up_xy*g_up_xy);

		auxfields.u_aux[SWAP_X]GRD_IND = 
		(g_up_xx*(*fields_ptr).u[M_C_X]GRD_IND + g_up_xy*(*fields_ptr).u[M_C_Y]GRD_IND)*Sdet;
		auxfields.u_aux[SWAP_Y]GRD_IND = 
		(g_up_yy*(*fields_ptr).u[M_C_Y]GRD_IND + g_up_xy*(*fields_ptr).u[M_C_X]GRD_IND)*Sdet;
     }
   }
 }

DERIV(grid_ptr,(struct field *)&auxfields.u_aux[SWAP_X],
         (struct field *)&dfields.du[D_X_M_C_X],&coorx);
DERIV(grid_ptr,(struct field *)&auxfields.u_aux[SWAP_Y],
         (struct field *)&dfields.du[D_Y_M_C_Y],&coory);


                /* derivatives of M_A_T */

DERIV(grid_ptr,(struct field *)&(*fields_ptr).u[M_A_T],
         (struct field *)&dfields.du[D_X_M_A_T],&coorx);

DERIV(grid_ptr,(struct field *)&(*fields_ptr).u[M_A_T],
         (struct field *)&dfields.du[D_Y_M_A_T],&coory);


DERIV(grid_ptr,(struct field *)&(*fields_ptr).u[M_B_T],
         (struct field *)&dfields.du[D_X_M_B_T],&coorx);

DERIV(grid_ptr,(struct field *)&(*fields_ptr).u[M_B_T],
         (struct field *)&dfields.du[D_Y_M_B_T],&coory);


DERIV(grid_ptr,(struct field *)&(*fields_ptr).u[M_C_T],
         (struct field *)&dfields.du[D_X_M_C_T],&coorx);

DERIV(grid_ptr,(struct field *)&(*fields_ptr).u[M_C_T],
         (struct field *)&dfields.du[D_Y_M_C_T],&coory);










#ifdef DISSIPATION


DISS(grid_ptr,(struct field *)&(*fields_ptr).u[M_A_T],
         (struct field *)&dfields.du[DISS_M_A_T]);
DISS(grid_ptr,(struct field *)&(*fields_ptr).u[M_A_X],
         (struct field *)&dfields.du[DISS_M_A_X]);
DISS(grid_ptr,(struct field *)&(*fields_ptr).u[M_A_Y],
         (struct field *)&dfields.du[DISS_M_A_Y]);

DISS(grid_ptr,(struct field *)&(*fields_ptr).u[M_B_T],
         (struct field *)&dfields.du[DISS_M_B_T]);
DISS(grid_ptr,(struct field *)&(*fields_ptr).u[M_B_X],
         (struct field *)&dfields.du[DISS_M_B_X]);
DISS(grid_ptr,(struct field *)&(*fields_ptr).u[M_B_Y],
         (struct field *)&dfields.du[DISS_M_B_Y]);

DISS(grid_ptr,(struct field *)&(*fields_ptr).u[M_C_T],
         (struct field *)&dfields.du[DISS_M_C_T]);
DISS(grid_ptr,(struct field *)&(*fields_ptr).u[M_C_X],
         (struct field *)&dfields.du[DISS_M_C_X]);
DISS(grid_ptr,(struct field *)&(*fields_ptr).u[M_C_Y],
         (struct field *)&dfields.du[DISS_M_C_Y]);
#endif

#ifdef DEBUG
   printf("in eqn. after taking derivatives \n");
#endif


/*************************************************************************************************/
/*************************************************************************************************/

/* inner points */

/*************************************************************************************************/

 {register int grid_ind1, grid_ind2;
 for (grid_ind1 = ni_1; grid_ind1 < nf_1; ++grid_ind1){
   for (grid_ind2 = ni_2; grid_ind2 < nf_2; ++grid_ind2){


#include"wave_2d_eq_eq.c"


	     (*derivs_ptr).u[M_A  ]GRD_IND =  m_a_t;
	     (*derivs_ptr).u[M_A_T]GRD_IND =  m_a_t_t;
	     (*derivs_ptr).u[M_A_X]GRD_IND =  m_a_x_t;
	     (*derivs_ptr).u[M_A_Y]GRD_IND =  m_a_y_t;

	     (*derivs_ptr).u[M_B  ]GRD_IND =  m_b_t;
	     (*derivs_ptr).u[M_B_T]GRD_IND =  m_b_t_t;
	     (*derivs_ptr).u[M_B_X]GRD_IND =  m_b_x_t;
	     (*derivs_ptr).u[M_B_Y]GRD_IND =  m_b_y_t;

	     (*derivs_ptr).u[M_C  ]GRD_IND =  m_c_t;
	     (*derivs_ptr).u[M_C_T]GRD_IND =  m_c_t_t;
	     (*derivs_ptr).u[M_C_X]GRD_IND =  m_c_x_t;
	     (*derivs_ptr).u[M_C_Y]GRD_IND =  m_c_y_t;



	 }
     }
  }
 
/* --------------------------------------------- MPI -------------------------------------- */
/* ---------------------------------- check fields have been received and sent ------------ */

#include"mpi_W.c"


/* ------------------------------------------------------------------------------------------------------*/
/* ------------------------------------SWAP order of some faces------------------------------------------*/
/* ------------------------------------------------------------------------------------------------------*/

#include"wave_2d_eq_inv.c"

/* ------------------------------------------------------------------------------------------------------*/
/* --------------------------------SWAP direction of some fields-----------------------------------------*/
/* ------------------------------------------------------------------------------------------------------*/
 
#include"wave_2d_eq_jaco.c"

/*-------------------------------------------------------------------------------------------------------*/
/*-------------------------------------------------------------------------------------------------------*/

 /*----------------------------- BOUNDARY CONDITIONS --------------------------*/

  /* Boundary conditions at x=xi and x=xf */

	 /************************* x = xi ******************************/

{register int grid_ind1;
 register int grid_ind2;
 
 grid_ind1 = ni_1;

       nx=nx_x0;
       ny=ny_x0;
 
  
 for (grid_ind2 = ni_2; grid_ind2 < nf_2; ++grid_ind2){

 m_a_t_n = sfield_recv.s_x_0_u[M_A_T][grid_ind2];
 m_a_x_n = sfield_recv.s_x_0_u[M_A_X][grid_ind2];
 m_a_y_n = sfield_recv.s_x_0_u[M_A_Y][grid_ind2];

 m_b_t_n = sfield_recv.s_x_0_u[M_B_T][grid_ind2];
 m_b_x_n = sfield_recv.s_x_0_u[M_B_X][grid_ind2];
 m_b_y_n = sfield_recv.s_x_0_u[M_B_Y][grid_ind2];

 m_c_t_n = sfield_recv.s_x_0_u[M_C_T][grid_ind2];
 m_c_x_n = sfield_recv.s_x_0_u[M_C_X][grid_ind2];
 m_c_y_n = sfield_recv.s_x_0_u[M_C_Y][grid_ind2];


#include"wave_2d_face.c"

     }
 }
   

	 /************************* x = xf ******************************/

{register int grid_ind1;
 register int grid_ind2;


      nx=nx_x1;
      ny=ny_x1;

      grid_ind1 = nf_1-1;

 for (grid_ind2 = ni_2; grid_ind2 < nf_2; ++grid_ind2){

 m_a_t_n = sfield_recv.s_x_1_u[M_A_T][grid_ind2];
 m_a_x_n = sfield_recv.s_x_1_u[M_A_X][grid_ind2];
 m_a_y_n = sfield_recv.s_x_1_u[M_A_Y][grid_ind2];

 m_b_t_n = sfield_recv.s_x_1_u[M_B_T][grid_ind2];
 m_b_x_n = sfield_recv.s_x_1_u[M_B_X][grid_ind2];
 m_b_y_n = sfield_recv.s_x_1_u[M_B_Y][grid_ind2];

 m_c_t_n = sfield_recv.s_x_1_u[M_C_T][grid_ind2];
 m_c_x_n = sfield_recv.s_x_1_u[M_C_X][grid_ind2];
 m_c_y_n = sfield_recv.s_x_1_u[M_C_Y][grid_ind2];


#include"wave_2d_face.c"





  }
  }
 


  /* Boundary conditions at y = yi and y = yf */

	 /************************* y = yi ******************************/

 {register int grid_ind1;
 register int grid_ind2;

       nx=nx_y0; 
       ny=ny_y0; 

       grid_ind2 = ni_2;

 for (grid_ind1 = ni_1; grid_ind1 < nf_1; ++grid_ind1){

 m_a_t_n = sfield_recv.s_y_0_u[M_A_T][grid_ind1];
 m_a_x_n = sfield_recv.s_y_0_u[M_A_X][grid_ind1];
 m_a_y_n = sfield_recv.s_y_0_u[M_A_Y][grid_ind1];

 m_b_t_n = sfield_recv.s_y_0_u[M_B_T][grid_ind1];
 m_b_x_n = sfield_recv.s_y_0_u[M_B_X][grid_ind1];
 m_b_y_n = sfield_recv.s_y_0_u[M_B_Y][grid_ind1];

 m_c_t_n = sfield_recv.s_y_0_u[M_C_T][grid_ind1];
 m_c_x_n = sfield_recv.s_y_0_u[M_C_X][grid_ind1];
 m_c_y_n = sfield_recv.s_y_0_u[M_C_Y][grid_ind1];




#include"wave_2d_face.c"

 }
 }


	 /************************* y = yf ******************************/

 {register int grid_ind1;
 register int grid_ind2;

      nx=nx_y1;
      ny=ny_y1;

	 grid_ind2 = nf_2-1;

 for (grid_ind1 = ni_1; grid_ind1 < nf_1; ++grid_ind1){

 
 m_a_t_n = sfield_recv.s_y_1_u[M_A_T][grid_ind1];
 m_a_x_n = sfield_recv.s_y_1_u[M_A_X][grid_ind1];
 m_a_y_n = sfield_recv.s_y_1_u[M_A_Y][grid_ind1];

 m_b_t_n = sfield_recv.s_y_1_u[M_B_T][grid_ind1];
 m_b_x_n = sfield_recv.s_y_1_u[M_B_X][grid_ind1];
 m_b_y_n = sfield_recv.s_y_1_u[M_B_Y][grid_ind1];

 m_c_t_n = sfield_recv.s_y_1_u[M_C_T][grid_ind1];
 m_c_x_n = sfield_recv.s_y_1_u[M_C_X][grid_ind1];
 m_c_y_n = sfield_recv.s_y_1_u[M_C_Y][grid_ind1];

#include"wave_2d_face.c"

 }
  }
 
  

 

/*  printf("At the end of eq. in n%d \n", grid); */
}
/*************************************************************************************/

 




